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1.  Introduction 


This  project  investigates  computational  modeling  of  fatigue  crack  growth  in  spiral 
bevel  gears.  Current  work  is  a  continuation  of  the  previous  efforts  made  to  use  the 
Boundary  Element  Method  (BEM)  to  simulate  tooth-bending  fatigue  failure  in  spiral 
bevel  gears.  This  report  summarizes  new  results  predicting  crack  trajectory  and  fatigue 
life  for  a  spiral  bevel  pinion  using  the  Finite  Element  Method  (FEM). 

Predicting  crack  trajectories  is  important  in  determining  the  failure  mode  of  a 
gear.  Cracks  propagating  through  the  rim  may  result  in  catastrophic  failure,  whereas  the 
gear  may  remain  intact  if  one  tooth  fails  and  this  may  allow  for  early  detection  of  failure. 
Being  able  to  predict  crack  trajectories  is  insightful  for  the  designer.  However,  predicting 
growth  of  three-dimensional  arbitrary  cracks  is  complicated  due  to  the  difficulty  of 
creating  three-dimensional  models,  the  computing  power  required,  and  absence  of  closed- 
form  solutions  of  the  problem. 

Previously  in  this  project,  tooth-bending  fatigue  failure  in  spiral  bevel  gears  was 
investigated  using  the  BEM  [1].  These  analyses  were  significant  in  developing  a  method 
for  predicting  three-dimensional,  non-proportional,  fatigue  crack  growth  incorporating 
moving  loads.  Prior  to  the  BEM  study,  there  were  not  many  examples  of  three- 
dimensional  crack  growth  work  on  gears.  Most  of  the  gear  studies  in  the  literature  are 
two-dimensional  analyses.  Furthermore,  moving  loads  had  not  been  considered  in  these 
analyses. 

The  BEM  predictions  were  helpful  to  determine  areas  in  need  of  further  research. 
These  were:  improving  accuracy  of  simulation  by  switching  from  the  boundary  element 
to  the  finite  element  method,  decreasing  computation  time  by  employing  parallel  FEM 
analysis,  and  improving  modeling  of  contact  tooth  loads. 

Substantial  improvement  in  computation  time  has  been  achieved  by  using  parallel 
FEM  technologies  on  PC-clusters.  This  makes  it  possible  to  simulate  crack  growth  in 
large  models,  like  the  current  gear  problem,  in  a  few  minutes.  It  is  also  possible  to  obtain 
accurate  values  of  stress  intensity  factors  (SIF)  using  the  FEM  and  elastic,  equivalent 
domain  J-integral  approaches.  A  new  technique  was  developed  for  modeling  moving 
tooth  loads,  which  decreases  the  computation  time  while  introducing  a  more  efficient 
procedure  for  applying  contact  loads. 

Another  focus  of  this  project  was  performing  three-dimensional  contact  analysis 
of  a  spiral  bevel  gear  set  incorporating  cracks.  These  analyses  were  significant  in 
determining  the  influence  of  change  of  tooth  flexibility  due  to  crack  growth  on  the 
magnitude  and  location  of  contact  loads.  This  is  an  important  concern  since  change  in 
contact  loads  might  lead  to  differences  in  SIFs  and  therefore  result  in  alteration  of  the 
crack  trajectory. 

Contact  analyses  performed  in  this  report  showed  the  expected  trend  of 
decreasing  tooth  loads  carried  by  the  cracked  tooth  with  increasing  crack  length. 
Decrease  in  tooth  loads  lead  to  differences  between  SIFs  extracted  from  finite  element 
contact  analysis  and  finite  element  analysis  with  Hertz  contact  loads.  This  effect  became 
more  pronounced  as  the  crack  grew. 
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2.  Three-dimensional  FEM  Simulations  of  Fatigue  Crack  Growth  in  Spiral  Bevel 
Gears 

The  main  objective  of  the  current  work  is  to  demonstrate  that  we  can  improve 
accuracy  and  efficiency  of  crack  growth  simulations  in  spiral  bevel  pinion  through  FEM 
computations  along  with  improved  modeling  of  tooth  loads.  In  the  current  study,  as  in 
previous  BEM  simulations,  LEFM  theories  are  used  for  fatigue  crack  growth  predictions 
in  gears.  Fatigue  crack  growth  rates  were  determined  using  a  modified  Paris  model 
accounting  for  crack  closure.  Crack  trajectory  predictions  are  made  using  the  maximum 
circumferential  stress  theory.  Moving  loads  were  taken  into  account  by  discretizing  the 
loading  into  15  steps  that  created  a  non-proportional  load  history  at  the  tooth  root. 

Accuracy  of  computations  can  be  improved  using  the  FEM  that  allows  a  better 
SIF  calculation  method,  the  equivalent  domain  J-integral  method.  Solution  time  for  each 
load  step  can  be  decreased  substantially  by  employing  parallel  FEM  analysis.  Another 
improvement  in  the  accuracy  of  the  results  can  be  obtained  by  developing  a  better 
representation  of  moving  contact  loads.  This  is  accomplished  by  approximating  the 
moving  contact  ellipse  by  interpolation  using  the  shape  functions  on  the  surfaces  of  the 
loaded  elements.  In  the  earlier  boundary  element  computations,  contact  ellipses  were  a 
part  of  the  geometry  models.  Since  the  ellipses  overlapped,  four  different  geometry 
models  were  needed  for  each  crack  configuration.  In  the  current  approach,  contact  loads 
are  no  longer  a  part  of  the  geometry  model.  This  approach  brings  a  substantial  decrease  in 
computation  time  since  only  one  geometry  model  for  each  crack  configuration  is  now 
needed. 

The  method  used  in  FEM  simulations  of  spiral  bevel  gear  is  composed  of  the 
following  steps: 

1 .  Create  initial  geometry  model  of  the  spiral  bevel  pinion  using  OSM. 

2.  Specify  boundary  conditions  and  material  properties  in  FRANC3D. 

3.  Initiate  crack  in  FRANC3D. 

4.  Create  a  surface  mesh  composed  of  triangular  elements  in  FRANC3D. 

5.  Create  a  3D  FEM  mesh  of  the  model  composed  of  tetrahedra  by  J-Mesh  using  a 
geometry  description  file  written  out  from  FRANC3D. 

6.  Calculate  the  magnitude  and  location  of  the  nodal  contact  loads  for  each  load  step 
on  the  loaded  elements  using  the  mesh  file  produced  by  J-Mesh. 

7.  Run  FEM  analysis  on  a  parallel  PC-cluster  for  15  discrete  load  steps.  Using  the 
displacement  results,  calculate  SIF  values  at  each  discrete  crack  front  point  for 
each  load  step. 

8.  Calculate  the  amount  of  crack  extension  at  each  discrete  crack  front  point  as  a 
result  of  15  steps  of  non-proportional  moving  load. 

9.  Determine  new  crack  front  by  piecewise  least  squares  fit  of  the  propagated  points 
corresponding  to  each  discrete  crack  front  point  in  FRANC3D. 

10.  Remesh  surface  locally  and  repeat  steps  5  to  10. 
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2.1  Modeling  of  Moving  Tooth  Loads 


Contact  in  a  spiral  bevel  gear  occurs  in  three-dimensions  following  a  path  along 
the  tooth  surface  starting  from  the  fillet  of  the  toe  to  the  top  of  the  heel.  In  analyzing  the 
spiral  bevel  gear,  it  is  assumed  that  the  contact  between  the  gear  and  pinion  follows  Hertz 
theory  of  elastic  contact.  Hertzian  contact  holds  as  long  as  the  significant  dimensions  of 
the  contact  area  are  small  compared  to  the  dimensions  of  each  body  and  to  the  relative 
radii  of  curvature  of  surfaces  [2].  Hertz  theory  assumes  that  surfaces  are  continuous  and 
non-conforming,  strains  are  small,  each  solid  can  be  considered  as  an  elastic  half-space 
and  surfaces  are  frictionless. 

Under  these  assumptions,  only  normal  pressure  acts  between  two  bodies  due  to 
Hertz  contact  producing  normal  displacements  of  the  surfaces.  Since  the  contact  area  is 
assumed  to  be  elliptical  the  region  is  defined  by, 


with  semi-ellipsoidal  pressure  distribution  of  the  form, 


f  2  2  V/2 

_ y_ 

P  Po  1  2  »2 

a  °  J 


Total  load  acting  on  the  ellipse  and  average  load  are  then  calculated  as, 

2nab  2 


Po 


avg 


Po 


(2) 

(3) 


In  the  above  equations  po  denotes  the  Hertz  stress  which  is  the  maximum  stress,  a  and  b 
are  half  of  major  and  minor  axes  of  the  ellipse,  respectively,  Figure  2.1. 


Figure  2.1:  Representation  of  Hertz  stress  over  an  elliptical  area. 

In  order  to  determine  the  contact  points  on  the  gear  tooth,  the  method  described  in 
[3]  is  used  to  map  the  three-dimensional  contact  path  into  a  two  dimensional  space.  The 
following  equations  were  used  for  the  transformation: 


PnMXI  ~  1,2  +  ^NMX2  ) 

4 n  =  ~  ZM )  cos  ex  +  (Rtf  ~  Rm  )  sin  ex 
r\N  =(Rn  -RM)cosa  +  (ZN  -ZM) sin  a 


sin  a  = 


^1  ^2 
2  a 


cos  a  = 


■^1  ^2 
2b 


(4) 
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Xm,Ym,Zm=  coordinates  of  mean  contact  point  (corresponds  to  the  contact  ellipse  centers) 
Xn,Yn,Zn  =  coordinates  of  an  arbitrary  contact  point 

Xu,Yi,2,Zu  =  coordinates  of  points  corresponding  to  1  and  2  in  Figure  2.2  in  three 
dimensions. 


k - - - w 

2a 


Figure  2.2:  Sketch  of  transformed  contact  ellipse  with  mean  contact  point  M  and  an 
arbitrary  contact  point  N. 

An  algorithm  was  developed  to  find  the  nodes  within  the  contact  ellipses  and  to 
calculate  the  applied  load  on  these  nodes.  This  algorithm: 

1 .  reads  the  nodal  information  from  the  mesh  and  geometry  files  created  by  3D  FEM 
meshing  program, 

2.  retrieves  the  geometry  information  of  the  n  th  contact  ellipse, 

3.  transforms  the  nodal  points  defined  in  3D  coordinates  into  contact  plane, 

4.  checks  if  the  mesh  nodes  are  within  the  contact  area  bounded  by  the  rt  th  ellipse 
(n  =  1 ...  1 5),  eg.  Figure  2.3, 

5.  calculates  nodal  contact  loads  if  the  node  is  within  the  contact  area,  and 

6.  writes  out  the  information  of  nodes  found  by  this  search  and  corresponding  loads. 


Figure  2.3:  Nodal  points  where  contact  loads  are  applied  for  load  step  11. 


2.2  Parallel  FEM 

Understanding  fracture  phenomena  in  materials  is  crucial  for  science  and 
engineering  applications  such  as  crack  propagation  in  spiral  bevel  gears.  However,  crack 
growth  simulations  are  complicated  and  computationally  expensive  in  realistic  structures. 
In  order  to  solve  realistic  3D  problems,  it  is  necessary  to  develop  fast  and  accurate 
computational  simulation  tools. 
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The  BEM  is  an  alternative  to  the  FEM  for  solving  large  problems.  The  BEM  is 
advantageous  since  it  only  requires  a  surface  mesh  and  has  fewer  equations  to  be  solved. 
However,  it  entails  a  folly  dense,  nonsymmetric  system  of  equations,  and  it  is  difficult  to 
extract  accurate  near-crack-front  stresses  and  displacements.  In  this  respect,  it  is  more 
desirable  to  use  FEM.  This  brings  the  challenge  of  implementing  fast  and  robust  parallel 
sparse  solvers. 

A  parallel  PC-cluster  in  Cornell  Theory  Center  was  used  for  the  development  of  a 
parallel  FEM  solver  for  crack  growth  simulations.  The  hardware  consists  of  32  nodes  of  a 
64  node  cluster  of  Dell  PowerEdge  1650  servers  with  a  Giganet  interconnect  running 
Microsoft  Windows  2000  Advanced  Server.  Each  node  contains  2  Intel  Pentium  III 
processors  clocked  at  1  GHz  and  2  GB  of  RAM.  Only  one  processor  per  node  is  used  for 
the  simulations. 

Parallel  FEM  simulation  of  crack  propagation  makes  it  possible  to  simulate  large 
problems  that  were  impossible  or  impractical  to  solve  previously.  The  objective  in 
developing  such  a  system  is  to  be  able  to  solve  a  real-life  problem  with  multiple  crack 
propagation  steps  with  1,000,000  degrees  of  freedom  in  an  hour.  Such  a  problem  with 
10,000  degrees  of  freedom  takes  about  100  hours  using  a  state-of-the-art  single  processor 
workstation  [4]. 

Crack  propagation  in  a  spiral  bevel  gear  stands  as  a  large,  realistic  engineering 
problem  that  incorporates  the  computational  difficulties  mentioned  above.  In  that  respect, 
the  development  of  parallel  FEM  simulations  is  a  significant  step  in  achieving  the  goal  of 
developing  a  practical  and  accurate  numerical  design  tool  for  gears. 

2.3  SIF  Calculations 

Complex  problems,  such  as  the  current  gear  problem,  require  computation  of 
mixed  mode  SIFs  due  to  the  combined  loading  that  they  receive.  In  the  literature,  several 
approaches  for  computing  SIFs  can  be  found.  Among  these  methods,  the  displacement 
correlation  method  [5-8]  and  equivalent  domain  J-integral  methods  [9-14]  are  widely 
used  in  numerical  computations. 

In  the  previous  BEM  simulations  of  the  spiral  bevel  gear,  the  displacement 
correlation  method  was  used  for  calculating  SIFs.  This  method  only  requires 
displacement  information  on  the  boundary  near  the  crack  front  and  is  computationally 
efficient.  Since  BEM  solves  for  displacement  information  only  on  the  boundaries,  this 
method  is  suitable  for  the  BEM  analysis.  The  displacement  correlation  method  makes  use 
of  numerical  solutions  for  displacements  in  combination  with  the  analytic  solutions  of 
SIFs.  This  method  is  only  as  accurate  as  the  computed  displacements  at  the  correlation 
points. 

The  equivalent  J-integral  method  is  adaptable  to  complex  problems  and  presents 
much  more  flexibility  in  the  type  of  applications  for  which  it  can  be  used.  We  used  the 
equivalent  domain  J-integral  approach  instead  of  displacement  correlation  in  the  FEM 
simulations  of  the  spiral  bevel  gear. 

J-integral  is  the  energy  release  expression  associated  with  crack  growth  defined  as 
an  integral  along  a  contour  that  encloses  a  crack  tip.  J-integral  is  defined  as  [15], 
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m 

J=\{Wnx-viJ-^-nj)dr 

where  W  is  the  strain  energy  density,  n,  is  the  unit  outward  normal  to  the  contour,  Oy  is 
the  traction  exerted  on  the  body  enclosed  by  F  and  crack  surface.  This  integral  is 
evaluated  along  an  arbitrary  contour  enclosing  the  crack  tip  and  is  path  independent. 


Figure  2.4:  Paths  of  J-integral  calculations  at  a  crack  tip. 

Numerical  evaluation  of  the  J-integral  involves  the  use  of  equivalent  domain 
integral  approach.  In  this  approach,  the  integration  along  the  contour  is  replaced  by  an 
integral  over  a  finite  size  domain.  The  expression  for  J-integral  then  becomes. 


J„  =  -a,  | “pnjqdA  k- IX  (6) 

r V  k 

A  continuous  function  q(x i,x2)  is  introduced  in  the  above  equation  in  order  to  replace  the 
integration  contour  T  by  Tj.  Function  q  has  a  unit  value  at  the  crack  tip  and  it  is  zero 
along  P,  Ti  and  T2.  The  variation  of  this  function  within  the  domain  is  arbitrary.  Using 
divergence  theorem  equation  6  takes  the  form, 


J,  =  j[W 


dq_ 

dxt 


-a 


dui  dq 
,J  dxk  dxj 


VA  -  J{ 


dW 

dx. 


d  fay  zr-])qdA 


dxj 


dxt 


(7) 


Note  that  the  second  term  of  this  integral  vanishes  for  elastic  problems.  SIFs  are 
calculated  within  the  vicinity  of  the  3D  crack  front  by  the  pointwise  values  of  J-integral. 
3D  cracks  can  be  represented  by  2D  plane  strain  fields  at  the  crack  front  points.  These 
fields  are  satisfied  only  within  close  proximity  of  the  crack  tip,  however  it  is  well  known 
that  finite  element  calculations  of  the  field  quantities  are  not  very  accurate  close  to  the 
crack  tip.  Equivalent  domain  integral  approach  circumvents  this  difficulty  by  replacing 
the  line  integral  by  a  finite  domain  integral.  The  relation  between  the  J-integral  and  SIFs 
is  given  for  an  isotropic,  linear  elastic  material  as  follows, 
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Jj=—  (K^+K/)  and  J2  =-2(^  +  1}  KK 

1  8G  8  G 

In  the  above  equations  J2  is  not  path  independent.  For  a  pure  mode  I  or  mode  II  problem, 
SIFs  can  be  calculated  directly  from  these  equations  however  under  mixed  mode  loading 
X/and  Ku  cannot  be  calculated  separately.  Bui  [16]  proposed  a  mixed  mode  separation 
method  by  decomposing  the  elastic  displacement  field  into  symmetrical  and  anti- 
symmetrical  parts.  This  results  in  separation  of  J  into  two  components  that  are  given  as 
follows  [16], 


1  —  V2  2 

J  =  Jl+J2,  Jj  =  — — — Kj2  , 


1  —  V2  2 

Jn 


i  ,  11 

U,  —  U,  +  M, 


(9) 


Examples  of  numerical  investigations  on  mixed  mode  problems  using  this  approach 
yielded  highly  accurate  results  [17,18]. 

Numerical  techniques  generally  provide  good  results  in  an  average  sense  over  the 
domain.  However,  governing  equations  and/or  boundary  conditions  are  only 
approximately  satisfied  at  any  one  point.  This  indicates  that  displacement  correlation 
SIFs  are  not  very  accurate  since  they  depend  on  the  displacement  values  at  a  point.  On 
the  other  hand,  J-integral  SIFs  give  better  results  due  to  considering  much  more  of  the 
results  domain  over  which  approximation  errors  tend  to  cancel. 

Two  models  with  known  analytic  solution  are  investigated  to  show  the  higher 
accuracy  of  the  FEM  calculations.  The  first  model  is  a  single,  internal,  flat,  penny¬ 
shaped  crack  in  an  infinite  body.  The  radius  of  the  penny-shaped  crack  is  0.1.  This  crack 
is  centered  in  a  block  with  dimensions  10x5x5. 


t  to 


Figure  2.5:  Flat,  internal  penny-shaped  crack  subjected  to  remote  tension. 

For  this  problem  Mode  II  and  Mode  III  SIFs  are  zero.  The  analytical  solution  for 
the  Mode  I  SIF  is  given  by  [19], 


(10) 
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Table  2.1  shows  the  comparison  of  solutions  from  the  BEM,  and  FEM  and  the  analytical 
solution.  BEM  results  shown  in  Table  2.1  are  calculated  by  the  displacement  correlation 
method  using  linear  quadrilateral  elements  with  12  subdivisions  along  the  radius  and 
along  the  quarter  of  the  crack  front.  FEM  calculations  are  carried  out  using  the  equivalent 
domain  J-integral  method  with  36,583  tetrahedra.  Note  that  an  elastic  modulus  of  10,000 
and  a  Poisson’s  ratio  of  0.0  are  assigned  to  the  model,  and  a  uniform,  far-field  tension  of 
unity  is  applied. 

Table  2. 1 :  Average  Mode  I  SIF  values  along  the  crack  front  calculated  by  the  BEM  and 
the  FEM  compared  with  the  analytical  solution. 


Solution  method 

Penny-shaped 

%  Error 

Analytical 

TTMl 

- 

BEM 

—.IBM 

2.6 

FEM 

0.352 

1.4 

The  second  example  is  a  single,  internal,  flat,  elliptical  crack  in  an  infinite  body 
with  an  elastic  modulus  of  10,000  and  a  Poisson’s  ratio  of  0.0.  The  orientation  of  the 
crack  is  the  same  as  in  Figure  2.5.  An  elliptical  crack  with  0.1  minor  and  0.2  major  axis 
dimensions  is  embedded  in  a  10  x  5  x  5  block.  A  uniform  far-field  tension  of  unity  is 
applied.  The  analytic  solution  of  Kj  involves  elliptic  integrals  whereas  A//  and  Km  are 
zero.  Detailed  discussion  of  the  analytical  solution  can  be  found  in  [19].  Comparison  of 
solutions  from  the  BEM,  and  FEM  and  the  analytical  solution  are  given  in  Table  2.2. 
BEM  calculations  are  done  using  the  displacement  correlation  method  with  linear 
quadrilateral  elements  with  5  subdivisions  along  the  axes  and  12  subdivisions  along  the 
quarter  of  the  crack  front.  The  FEM  results  are  computed  using  J-integral  method  with 
13,464  tetrahedra. 

Table  2.2:  Maximum  and  minimum  Mode  I  SIF  values  along  the  crack  front  calculated 
by  the  BEM  and  the  FEM  compared  with  the  analytical  solution. 


Solution  method 

Max  A} 

%  Error 

Min  A/ 

%  Error 

Analytical 

0.463 

- 

- 

BEM 

0.435 

5.9 

0.304 

7.0 

FEM 

0.442 

4.5 

1.8 

2.4  Determining  New  Crack  Front 

Crack  extension  due  to  one  non-proportional  load  cycle  is  computed  using  the 
approach  developed  in  [1].  As  previously  mentioned,  a  modified  Paris  model  was  used 
that  accounts  for  crack  closure  effects.  A  crack  is  assumed  to  advance  when  its  SIF  is 
large  enough  to  overcome  closure  and  is  larger  than  the  SIF  of  the  previous  load  step. 
The  new  crack  front  is  determined  as  follows: 
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1 .  For  every  point  along  the  crack  front: 

a)  Calculate  amount  of  crack  growth,  da'  corresponding  to  each  load  step 
1, — 15. 

K/-K/-1 


da '  =• 


A/C 


da 


da 

dN 


=  C(AK^y 


A KeJf=Kr-K, 


op 


K. 


op 


K„ 


=  A0  +  A^R  +  AjR*  +  AgR3 


Ag  =(0.825-  0.34a  +  0.05a2) 


4  =  (0.415a-  0.071a2 


cos(^) 

2an 


— il/cx 


for  R>0 


A2  —  1  Ag  A^  Ag 
Ag  =  2  Ag  +  A j  —  1 


2. 

3. 


b)  Calculate  the  angle  of  crack  growth  corresponding  to  each  load  step  using 
maximum  circumferential  stress  theory 


0'  =  2  tan 


i*/±i 

4  Ku‘  4^| 


(*4)-+. 

K  ' 


II 


c)  Calculate  the  final  coordinates  of  the  crack  front  and  trajectory  angles  by 
approximating  the  contributions  from  each  load  step  by  a  straight  line. 

d)  Determine  number  of  cycles,  N,  necessary  to  have  a  reasonable  amount  of 
crack  advance  compared  to  pinion’s  geometry. 

Update  geometry  using  the  new  crack  front  points. 

Repeat  1  and  2  for  each  crack  growth  step. 


2.5  Results 

Three-dimensional  finite  element  analyses  were  performed  for  an  OH-58  spiral 
bevel  pinion  for  54  crack  growth  steps,  Figure  2.6.  In  the  previous  BEM-based  analysis, 
only  13  growth  steps  were  performed.  One  crack  growth  step  consists  of  15  static  finite 
element  analyses  in  order  to  simulate  moving  loads  on  the  pinion  tooth.  After  the  54 
crack  growth  step  the  analysis  was  stopped  because  it  was  concluded  that  propagating  the 
crack  further  would  not  provide  additional  insights  into  the  simulations  and  results.  All 
predictions  were  based  on  the  methodology  outlined  in  section  2. 
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Figure  2.6:  A  spiral  bevel  pinion  model  with  surface  FE  mesh  (a)  whole  pinion  (b)  close- 
up  view  of  pinion  teeth.  Note  the  location  of  the  initial  crack  in  the  middle  tooth. 


Figures  2.7  through  Figure  2.9  show  Mode  I,  II  and  III  stress  intensity  factors  for 
the  initial  crack  configuration  for  the  first  eleven  load  steps.  Load  steps  1-4  correspond  to 
double  tooth  contact,  load  steps  5-11  correspond  to  single  tooth  contact  and  load  steps 
12-15  again  correspond  to  double  tooth  contact.  Since  a  crack  is  not  assumed  to  advance 
when  its  SIF  is  smaller  than  the  SIF  of  the  previous  load  step,  only  the  first  11  steps 
practically  contribute  to  crack  growth  calculations.  In  these  figures  crack  front  position 
corresponds  to  the  points  near  the  crack  front  at  which  SIF  values  are  evaluated.  In  the 
initial  step  of  crack  growth  simulations  the  crack  front  was  discretized  such  that  there 
were  222  points  at  which  SIF  values  were  calculated.  This  is  compared  to  only  61  points 
previously  used  with  the  BEM  model  [1]. 

Figure  2.10  shows  the  crack  trajectory  predictions  by  the  FEM  on  the  tooth 
surface  and  cross  section  of  the  tooth  for  several  crack  growth  steps  including  the  initial 
and  final  configurations  of  the  crack.  Figure  2.11  is  a  close-up  view  of  the  last  step  of 
crack  growth  from  the  toe  end  of  the  tooth. 


—  Load  step  1  Load  step  3  —  Load  step  5  —  Load  step  7  —  Load  step  9  — -  Load  step  1 1  j 

Figure  2.7:  Initial  Mode  I  SIF  distribution  for  load  steps  one  through  eleven.  Note  that 
orientation  is  from  heel  to  toe. 
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Figure  2.8:  Initial  Mode  II  SIF  distribution  for  load  steps  one  through  eleven.  Note  that 
orientation  is  from  heel  to  toe. 


Figure  2.9:  Initial  Mode  HI  SIF  distribution  for  load  steps  one  through  eleven.  Note  that 
orientation  is  from  heel  to  toe. 
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Initial  Step:  Step  0: 


Step30: 


/"\ 


Figure  2.10:  Crack  trajectory  predictions  by  FEM  on  tooth  surface  and  tooth  cross  section. 
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Figure  2.11:  Close-up  view  of  the  crack  after  the  last  step  of  crack  growth. 

2.6  Comparison  of  FEM  and  BEM  Results 

This  section  compares  FEM  and  BEM  results  in  terms  of  accuracy  and  efficiency. 
BEM  results  presented  in  this  section  are  from  reference  [1].  One  of  the  most  significant 
improvements  in  the  FEM  approach  over  the  BEM  simulations  is  the  difference  in  the 
solution  time  of  the  problem.  This  comparison  is  summarized  in  Table  2.3. 

Table  2.3:  Comparison  of  performance  between  FEM  and  BEM  simulations. 


Results 

FEM 

BEM 

Computer  power 

32  nodes* 
PC-cluster,  1  GHz 

Single  DEC-Alpha 
Workstation,  175  MHz 

Number  of  unknowns 

-1-3  million 

6700 

(coefficient  matrix) 

(sparse,  symmetric) 

(dense,  unsymmetric) 

Solution  time  per  load  cycle 

5-15  hours 

~60  hours 

(1  load  cycle=15  load  steps) 

(20min-l  h/load  step) 

(4  models) 

Typical  time  for  one  crack  step 

~8-18  hours* 

—65  hours* 

Solution  time  plus  3  hours  post-processing  and  remeshing 
Solution  time  plus  5  hours  post-processing  and  remeshing 


The  performance  comparison  in  Table  2.3  clearly  shows  that  there  is  a  very 
significant  improvement  in  computation  time  with  the  development  of  a  parallel  FEM 
solver.  The  solution  time  for  BEM  simulations  was  about  65  hours  for  one  crack  growth 
step  whereas  FEM  simulations  took  much  less  time  for  each  crack  growth  step  even 
though  the  models  have  many  more  degrees  of  freedom.  The  number  of  unknowns  for  the 
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FEM  was  about  a  million  at  the  initial  step  and  increased  to  more  than  3  million  as  the 
crack  front  advanced.  3D  volume  meshing  of  one  model  ranged  from  20  minutes  for  the 
smallest  model  to  about  an  hour  for  the  largest  models. 

Improvement  in  the  computation  time  for  each  crack  growth  step  allowed  us  to 
increase  solution  resolution  by  taking  smaller  crack  increment  steps  and  to  simulate  more 
steps  of  propagation.  Figures  2.12  to  2.14  show  the  comparison  of  crack  area,  depth  and 
crack  front  length  by  FEM  and  BEM.  These  plots  indicate  that  the  FEM-predicted  fatigue 
growth  rate  is  marginally  lower  when  we  compare  the  crack  area  and  the  crack  front 
length  and  more  cycles  are  required  to  achieve  the  same  amount  of  crack  growth  when 
compared  with  previous  BEM-based  results.  Figures  2.15  to  2.17  show  the  SIF 
comparisons  between  FEM  and  BEM  for  load  steps  3,  7,  and  11. 

Table  2.4:  Comparison  of  crack  growth  amount  between  FEM  and  BEM  simulations. 


Results 

FEM 

BEM 

Number  of  steps 

54 

13 

Total  number  of  cycles 

383,000 

311,000 

Final  crack  front  length 

1.55  in 

1.42  in 

Final  crack  area 

0.319  in2 

0.186  in2 

Final  crack  depth 

0.268  in 

0.188  in 

Number  of  cycles 

Figure  2.12:  Comparison  between  crack  area  versus  number  of  cycles  predicted  by  BEM 
and  FEM  analyses. 
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Figure  2.13:  Comparison  between  crack  depth  versus  number  of  cycles  predicted  by 
BEM  and  FEM  analyses. 
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18000 


Load  Step  3 


13000  H 


-7000 


Crack  front  position  (in) 

Mode  I,FEM  •-<•••  Mode  I,BEM  — Modell.FEM 
-  -  o  -  -  Moden,BEM  — <>—  Mode  III,FEM  -  -  o  Mode  m3EM 


Figure  2.15:  Comparison  of  initial  Mode  I,  II  and  III  SIFs  computed  by  the  FEM  and  the 
BEM  for  load  step  3. 


— *—  Mode  I,FEM  Mode  I3EM  — 13 —  Mode  II,FEM 

-  o  -  •  Mode  II,BEM  — ° —  Mode  III,FEM  --<>■■  Mode  III,BEM 

Figure  2.16:  Comparison  of  initial  Mode  I,  II  and  III  SIFs  computed  by  the  FEM  and  the 
BEM  for  load  step  7. 
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18000 


Load  step  1 1 


Crack  front  position  (in) 


— Mode  I,FEM  Mode  I,BEM  — o—  Mode  IIJEM 

•  •  o  -  -  Mode  n,BEM  — e-  Mode  III, FEM  -  -<>  -  •  Mode  m,BEM 

Figure  2.17:  Comparison  of  initial  Mode  I,  II  and  IE  SEFs  computed  by  the  FEM  and  the 
BEM  for  load  step  1 1. 


—♦-Mode  I  JEM  — *  -  -  Mode  I,BEM  — Mode  II.FEM 
-o -  Mode II JEM  Mode III,FEM  Mode  IIIJEM 


Figure  2.18:  Comparison  of  typical  SIF  history  computed  by  the  FEM  and  the  BEM  for  a 
load  cycle  at  the  midpoint  of  the  initial  crack  front. 
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Figure  2.19:  Comparison  of  typical  non-proportional  load  cycle  at  the  midpoint  of  the 
initial  crack  front  computed  by  the  FEM  and  the  BEM. 

Figures  2.18  and  2.19  present  the  comparison  of  SIF  variation  over  one  load  cycle 
and  non-proportional  load  history.  In  Figure  2.18,  it  is  seen  that  Kj  predictions  are  lower 
in  the  FEM  analysis  than  the  BEM  analysis  over  most  of  the  steps.  As  seen  in  Figure 
2.19,  the  ratio  of  Kn  to  Kj  over  15  load  steps  shows  a  similar  trend  as  observed  in  the 
BEM  analysis  with  a  slight  difference  in  values. 


Figure  2.20:  Comparison  of  crack  trajectory  predicted  by  the  BEM,  the  FEM  and  an 
experiment  on  the  tooth  surface  and  tooth  cross  section  at  the  midpoint  of  crack. 
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Crack  trajectory  comparisons  in  Figure  2.20  show  that,  on  the  tooth  surface,  crack 
trajectory  predicted  by  FEM  is  similar  to  the  BEM  results.  The  deviation  of  crack  path  at 
the  toe  end  in  the  BEM  predictions  was  also  observed  in  FEM  predictions.  On  the  heel 
end,  FEM  results  exhibit  a  less  steep  propagation  angle  than  the  experimental  trajectory. 
On  the  tooth  cross-section  it  is  observed  that  the  predicted  crack  trajectory  by  FEM  is  less 
steep  than  BEM  and  experimental  results.  FEM  results  also  exhibit  a  shallower  ridge 
formation  than  the  experiment.  Although  the  FEM-predicted  crack  path  still  does  not 
exactly  represent  the  crack  trajectory  observed  in  experiments,  it  captures  the  global 
behavior  of  the  crack  propagation  trajectory. 


3.  Three-dimensional  Finite  Element  Contact  Analysis 

The  meshing  of  a  gear  and  pinion  produces  loads  on  the  tooth  surface  that  change 
in  magnitude  and  position  in  time.  In  section  2.1  analyses  assuming  Hertzian  contact 
theory  were  described.  This  assumption,  however,  is  a  simplification  for  a  number  of 
reasons:  1)  the  teeth  and  rims  that  transmit  the  load  are  flexible,  2)  assumed  elliptical 
contact  area  on  the  tooth  surface  may  not  be  accurate  since  the  curvatures  of  the  teeth 
surfaces  in  contact  are  not  constant  over  the  contact  zone,  3)  the  center  of  contact,  where 
the  maximum  contact  stress  occurs,  may  differ  from  the  theoretical  contact  point,  and  4) 
most  importantly  for  the  current  investigation,  the  tooth  flexibility  changes  as  the  crack 
propagates  in  a  tooth.  The  change  in  flexibility  due  to  cracks  may  shift  the  contact  area, 
and  change  the  magnitude  of  the  contact  loading. 

In  this  section  a  series  of  three-dimensional  finite  element  analyses  that  model  the 
contact  conditions  between  the  gear  and  pinion  explicitly  are  described.  The  purpose  of 
this  study  is  to  assess  what  influence  the  details  of  the  load  transfer  through  the  contact 
conditions  have  on  the  computed  stress  intensity  factors,  and  through  them  the  predicted 
crack  trajectory. 

Three-dimensional  finite  element  contact  analyses  in  spiral  bevel  gears  have  been 
undertaken  by  several  researchers  to  study  the  stresses  induced  at  the  root  and  on  the 
surface  of  the  tooth.  Earlier  work  by  Bibel  et  al.  [20]  utilized  gap  elements  to  simulate  the 
contact  boundary  conditions  and  to  evaluate  the  stresses  on  the  tooth  surface.  Further 
study  on  this  subject  eliminated  the  need  of  using  gap  elements  through  the  use  of  contact 
capabilities  of  commercial  finite  element  programs.  Preliminary  results  showing  the 
moving  contact  patches  over  the  tooth  surface  are  reported  in  Bibel  and  Handschuh  [21]. 
These  analyses  also  incorporate  automatic  meshing  of  the  gear  and  pinion  via  the 
addition  of  user  subroutines  to  a  commercial  finite  element  program.  The  same  authors 
performed  additional  analyses  in  order  to  compare  experimentally  measured  tooth 
bending  stresses  of  spiral  bevel  gears  to  FEM  results  [22].  Other  studies  on  the  contact 
analysis  of  spiral  bevel  gears,  including  the  development  of  a  methodology  combining  a 
surface  integral  and  a  finite  element  solution,  are  described  by  Vijakaar  [23].  Numerical 
examples  carried  out  using  this  approach  show  the  change  in  the  contact  area  as  a  result 
of  the  misalignment  of  the  gears. 
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3.1  Approach 


The  goal  of  the  present  study  is  to  perform  contact  analyses  incorporating  cracks 
to  determine  the  effect  of  changing  tooth  flexibility  on  the  magnitude  and  location  of  the 
contact  load  and  therefore  on  the  stress  intensity  factors.  A  commercial  finite  element 
program,  ABAQUS  [24],  was  used  to  perform  the  contact  analysis  in  conjunction  with  the 
software  developed  by  the  CFG  to  calculate  fracture-related  parameters.  A  general 
outline  of  the  approach  is  as  follows: 

1 .  Create  initial  geometry  model  of  spiral  bevel  pinion  and  gear  with  OSM  using 
the  geometry  information  provided  by  NASA  Glenn  Research  Center. 

2.  Specify  boundary  conditions,  material  properties,  and  contact  surfaces  in 
FRANC3D. 

3.  Create  a  surface  mesh  composed  of  triangular  elements  in  FRANC3D. 

4.  Create  a  3D  FE  mesh  of  the  model  composed  of  tetrahedra  using  the  J-mesh 
program. 

5.  Transform  the  FRANC3D  model  into  an  ABAQUS  input  file. 

6.  Perform  contact  analysis  in  ABAQUS  and  determine  the  location  and  magnitude 
of  the  contact  loads. 

7.  Calculate  SlFs  using  the  results  of  FEM  contact  analysis  and  also  using  the  results 
from  analysis  with  Hertzian  contact  loads.  Compare  these  values  and  investigate 
their  effects  on  the  crack  trajectory  predictions. 

Geometries  of  a  one-tooth  sector  of  the  spiral  bevel  pinion  and  gear  were  created 
separately  in  OSM  using  the  information  provided  by  NASA.  This  information  regarding 
the  geometric  positions  of  points  defining  the  gear  set  was  determined  analytically  using 
the  method  developed  by  Litvin  et  al.  [3].  The  one-tooth  gear  and  pinion  models  were 
created  in  a  common  coordinate  system  with  the  same  axis  of  rotation  (z-axis). 
Additional  teeth  of  the  gear  and  pinion  were  created  by  copying  and  rotating  the  original 
tooth  by  an  angle  equal  to  360°  divided  by  number  of  gear  teeth.  The  resulting  models 
were  further  rotated  to  orient  them  in  the  proper  meshing  position.  At  the  meshing 
orientation  there  was  no  interference  between  the  gear  and  pinion  teeth. 

The  contact  analysis  was  performed  using  the  contact  pair  approach  in  ABAQUS. 
Contact  pairs,  named  as  master  and  slave  surfaces,  define  surfaces  which  can  potentially 
come  into  contact.  If  one  of  the  contact  surfaces  is  stiffer  or  has  a  coarser  mesh,  that 
should  be  selected  as  the  master  surface.  In  this  formulation  contact  direction  is  always 
normal  to  the  master  surface.  In  ABAQUS,  contact  can  occur  in  the  form  of  small-sliding 
or  finite-sliding.  Analyses  performed  in  this  work  used  the  finite-sliding  formulation 
where  surfaces  are  allowed  to  undergo  finite  separation  and  sliding.  This  formulation 
defines  the  master  surface  as  composed  of  element  faces.  Continuity  of  the  surface 
normals  was  attained  by  smoothing  the  transitions  between  the  elements.  Contact  surface 
tractions  were  calculated  assuming  that  the  surfaces  are  perfectly  hard,  which  means 
there  will  be  no  penetration  of  surfaces  into  one  another.  Contact  tractions  were  defined 
by  the  local  basis  system  formed  by  the  normal  to  the  master  surface. 

The  analyses  used  both  FRANC3D  and  ABAQUS.  FRANC3D  was  used  to  generate 
the  complicated  geometry  of  the  gear,  to  specify  the  analysis  attributes,  and  to  define  the 
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crack.  Also,  it  was  used  as  a  post-processor  to  calculate  the  SIFs.  ABAQUS  was  utilized 
for  finite  element  contact  analysis  of  the  model.  In  order  to  perform  contact  analysis  on 
models  with  cracks,  features  that  complement  FRANC3D  and  ABAQUS  were  needed.  In 
this  respect  FRANC3D  was  extended  to  include  a  capability  of  specifying  the  master  and 
slave  contact  surfaces  and  to  include  a  translator  to  convert  the  information  regarding  the 
geometry,  mesh,  boundary  conditions,  material  properties  and  contact  surfaces  to  an 
ABAQUS  input  file.  Also,  equivalent  domain  J-integral  implementation  was  added  to 
FRANC3D  to  obtain  more  accurate  SIF  values. 

Contact  analyses  involve  nonlinear  geometric  analysis.  The  gear/pinion  model  is 
meshed  with  linear  tetrahedra.  Volume  meshing  is  done  using  J-mesh  as  before.  The 
mesh  is  finer  on  the  contacting  tooth  surfaces  of  both  gear  and  pinion  compared  to  the 
rest  of  the  model.  In  performing  the  analyses,  it  was  assumed  that  there  was  no  friction 
between  the  contacting  surfaces. 

3.2  Analysis 

Finite  element  contact  analysis  of  a  spiral  bevel  gear  and  pinion  was  performed  on 
a  four-tooth  model  composed  of  two  pinion  teeth  and  two  gear  teeth.  This  model  is 
shown  in  Figure  3.1  and  is  adapted  from  [22].  Only  four  teeth  of  the  gear  and  pinion  were 
explicitly  modeled  since  this  is  the  minimum  number  of  teeth  that  is  needed  to  simulate 
double  tooth  contact.  We  chose  to  model  only  a  small  portion  of  the  gear  and  pinion  to 
decrease  the  computation  time  and  effort.  This  model  was  adequate  to  simulate  the 
contact  conditions  and  also  demonstrate  how  the  stress  intensity  factors  change  with 
increasing  crack  length. 

The  model  was  designed  such  that  load  applied  at  a  known  distance  from  the  axis 
of  rotation  transmits  the  required  amount  of  torque.  The  gear  model  was  extended  to 
include  a  solid,  "wedge  shaped"  part.  This  portion  was  a  modeling  artifact  added  so  that 
the  gear  could  be  allowed  to  rotate  freely  about  its  axis  of  rotation  and  come  into  contact 
with  the  pinion  as  a  result  of  the  applied  torque.  As  seen  in  Figure  3.1a,  the  pinion  teeth 
were  fixed  with  zero  displacements  in  x,y,z  directions  on  the  rim  surfaces.  Displacement 
components  of  the  gear  were  fixed  in  x,y,z  directions  along  the  line  of  axis  of  rotation. 
Gear  rim  surfaces  were  only  constrained  in  the  z  direction;  x  and  y  degrees  of  freedom  of 
the  gear  were  not  fixed  in  order  to  allow  for  free  rotation  in  that  plane.  Dimensions  of  the 
gear  sector  were  chosen  such  that  the  bending  effects  caused  by  the  application  of  load  is 
minimized.  Early  studies  showed  that  when  smaller  portions  of  the  gear  were  modeled 
contact  loads  changed  with  the  change  in  the  location  of  the  applied  load  due  to  the 
excessive  deformation  occurring  in  the  gear  hub. 

Full  design  torque  for  the  pinion  was  3099  in-lb.  In  order  to  obtain  the  torque  that 
is  transferred  by  the  gear  this  value  was  multiplied  by  the  gear  ratio,  giving  11580  in-lb. 
A  point  load  was  applied  at  a  distance  3.58  in.  from  the  gear  axis  of  rotation  to  produce 
this  amount  of  torque.  The  location  of  the  applied  load  is  shown  in  Figure  3.1a. 
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Figure  3.1:  Contact  model  of  spiral  bevel  gear  set.  (a)  whole  model  (b)  close  up  view  of 
gear  and  pinion  teeth  in  mesh  (c)  FE  mesh  of  the  model.  Note  that  two  teeth  of  pinion 
and  gear  are  explicitly  modeled.  Note  also  the  refined  mesh  on  the  pinion  and  gear  teeth 
surfaces  where  contact  is  expected  to  occur. 
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The  original  configuration  of  the  gear  and  pinion  corresponds  to  a  double  tooth 
contact  load  transfer.  In  order  to  obtain  a  single  tooth  contact  configuration  the  gear  set 
was  rotated  equivalent  to  a  12°  pinion  rotation.  During  the  analyses,  first  a  very  low 
displacement  was  applied  to  the  gear  in  order  to  ensure  contact  between  gear  and  pinion 
and  to  fix  the  model  in  space.  Then  the  load  was  increased  to  its  full  design  value. 

In  order  to  obtain  reasonable  contact  loads,  an  adequately  refined  mesh  was  very 
important.  Since  the  gear  has  a  doubly  curved  geometry,  a  fine  discretization  was  needed 
to  model  the  curved  edges.  In  our  studies,  we  observed  that  coarse  meshes  lead  to 
inaccurate  contact  stresses  and  stress  concentrations.  To  this  end,  we  discretized  the 
surfaces  where  contact  is  expected  to  occur  fine  enough  to  obtain  reasonable  results 
without  excessive  computational  cost. 

3.3  Results 

Contact  analyses  were  carried  out  on  four  models  which  have:  no  crack,  an  initial 
crack,  and  two  crack  sizes,  crack  #1  and  crack  #2.  Crack  #1  is  within  the  fillet  region  of 
the  pinion  tooth.  Crack  #2  is  a  larger  crack  extending  to  the  surface  of  the  tooth.  These 
four  geometries  were  analyzed  for  double  tooth  contact  and  single  tooth  contact 
configurations  that  correspond  to  0°  and  12°  pinion  rotation,  respectively. 

The  number  of  nodes  and  elements  for  the  analyzed  models  ranged  between 
12,000-16,000  and  52,000-70,000,  respectively.  Solution  time  of  a  finite  element  contact 
analysis  is  about  30  minutes  on  a  1  GHz  Intel  Pentium  El  machine. 


(c) 

Figure  3.2:  Pinion  teeth  with  a)  initial  crack  b)  crack  #1  c)  crack  #2.  Note  that  in  all  the 
models,  crack  is  in  tooth  1. 
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3.3.1  Comparison  of  contact  load  and  area  calculated  by  Hertz  theory  and  three- 
dimensional  FEM  contact  analysis 

Results  of  contact  analyses  and  how  they  compare  to  the  Hertz  theory  are 
summarized  in  Table  3.1.  Total  contact  load  calculated  by  finite  element  analysis 
compared  well  with  the  total  load  predicted  by  Hertz  theory.  However,  distribution  of  this 
total  contact  load  between  two  adjacent  teeth  shows  some  discrepancies.  There  is  nearly  a 
10%  difference  in  the  values  between  contact  analysis  and  Hertz  theory  predictions  on 
each  tooth. 

We  were  unable  to  attain  single  tooth  contact  conditions  using  this  model.  At  all 
rotations,  there  is  always  at  least  a  very  little  portion  carried  by  the  adjacent  tooth.  This 
might  be  due  to  the  inaccuracies  introduced  by  the  geometric  model  of  the  gears.  The 
closest  case  to  single  tooth  loading  was  obtained  when  a  12°  pinion  rotation  was 
introduced  to  the  model.  This  case  was  assumed  to  correspond  to  load  step  1 1  during 
which  the  highest  contact  load  occurs. 

Contact  loads  calculated  by  finite  element  analyses  decrease  on  tooth  1,  where  the 
crack  is  located,  and  increase  on  tooth  2  as  the  crack  grows.  These  results  indicate  that  as 
the  crack  becomes  larger,  the  flexibility  of  the  pinion  tooth  is  enhanced,  therefore  the 
contact  loads  taken  up  by  the  two  adjacent  teeth  are  redistributed.  A  larger  portion  of  the 
total  contact  load  has  to  be  carried  by  the  uncracked  neighboring  tooth.  This  trend  is  in 
agreement  with  the  expected  behavior. 


Table  3.1:  Comparison  of  contact  loads  and  areas  between  FEM  contact  simulations  and 
Hertzian  stress  calculations  for  two  contact  configurations  with  different  crack  lengths. 
Note  that  in  the  first  configuration  the  Hertz  contact  load  given  below  corresponds  to  load 
step  3  for  tooth  1  and  load  step  13  for  tooth  2,  and  in  the  second  configuration  it 
corresponds  to  load  step  11. 


Model 

Contact  Load  (lb) 
Tooth  #1  Tooth  #2 

Contact  Area  (in  ) 
Tooth  #1  Tooth  #2 

Configuration  # 

1  ( Rotation  0 °) 

0.0126 

No  crack 

1800 

1399 

0.0129 

Initial  crack 

1794 

1404 

0.0129 

0.0126 

Crack  #1 

1792 

1406 

0.0123 

0.0126 

Crack  #2 

1743 

1457 

0.0177 

0.0126 

Hertz  Solution 

1584 

1602 

0.0114 

0.0146 

Configuration  # 

2  ( Rotation  12°) 

0.0018 

No  crack 

2899 

318 

0.0193 

Initial  crack 

2885 

322 

0.0193 

0.0018 

Crack  #1 

2881 

326 

0.0193 

0.0018 

Crack  #2 

2860 

348 

0.0252 

0.0018 

Hertz  Solution 

3196 

- 

0.0216 

“ 
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Figure  3.3:  Contact  pressures  in  on  pinion  tooth  surface  for  a  double  tooth  contact 
configuration  with  (a)  initial  crack  (b)  crack  #1  (c)  crack  #2.  Note  that  the  values  on  the 
contour  bars  are  in  psi. 


For  earlier  steps  of  the  crack  growth,  when  the  crack  is  still  in  the  fillet  of  the 
pinion  tooth,  as  in  crack  #1,  there  is  no  significant  change  in  the  contact  load  and  area  on 
the  tooth  as  seen  in  Table  3.1.  However,  when  the  crack  extends  into  the  tooth  surface,  as 
in  crack  #2,  its  effect  on  the  flexibility  of  the  tooth  is  more  pronounced.  This  results  in  a 
more  significant  decrease  in  the  contact  load  on  tooth  1 . 

There  is  no  significant  difference  in  the  centroid  of  contact  loads  calculated  by 
finite  element  contact  analysis  for  different  crack  lengths.  However,  locations  of  contact 
loads  determined  by  these  analyses  were  shifted  towards  the  toe  compared  to  the  Hertz 
theory  predictions. 
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Figure  3.4:  Contact  pressures  on  pinion  tooth  surface  for  an  almost  single  tooth  contact 
configuration  with  (a)  initial  crack  (b)  crack  #1  (c)  crack  #2.  Note  that  the  values  on  the 
contour  bars  are  in  psi. 
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Contact  areas  shown  in  Figures  3.3  and  3.4  are  almost  elliptical,  showing  the 
same  trend  as  the  Hertz  contact  predictions.  Also,  the  distribution  of  the  contact  loads  is 
close  to  an  ellipsoid  shape  predicted  by  Hertzian  contact.  Maximum  contact  loads 
calculated  by  contact  analysis  were  higher  than  the  predictions  by  the  Hertz  theory. 

Results  of  contact  analyses  are  very  sensitive  to  the  mesh  size  of  the  contacting 
surfaces.  If  a  coarse  mesh  is  used,  computed  contact  pressure  is  not  continuous  on  the 
surface  since  only  certain  nodes  are  in  contact.  In  order  to  overcome  this  issue,  the  mesh 
on  tooth  surfaces  was  refined  several  times  until  a  continuous  contact  pressure 
distribution  was  obtained.  The  slightly  wavy  pattern  observed  in  the  contact  areas  in 
Figures  3.3  and  3.4  are  an  artifact  of  the  mesh  size,  shape  and  orientation.  Meshes  on  the 
gear  and  pinion  tooth  surfaces  are  built  such  that  there  is  compatibility  of  the  mesh  nodes 
between  these  two  contact  surfaces.  This  is  important  since  the  contact  is  calculated 
based  on  the  interaction  of  mesh  nodes  between  contacting  surfaces.  If  the  nodes  of  one 
of  the  surfaces  fall  between  the  two  nodes  defining  an  edge  of  the  other  surface,  then  the 
contact  load  is  taken  up  by  either  of  the  nodes  that  might  lead  to  different  contact 
patterns.  The  discontinuous  contact  load  pattern  observed  on  crack  #2  (Figure  3.3c  and 
3.4c)  for  both  configurations  might  be  due  to  the  surface  mesh  of  the  tooth.  At  this  crack 
step,  the  crack  was  extended  to  the  surface  of  the  pinion  tooth.  This  distorted  the 
regularity  of  the  mesh  on  the  pinion  surface.  Since  the  compatibility  of  the  meshes  on 
both  surfaces  is  very  important  for  the  contact  analysis,  the  loss  of  this  compatibility 
might  lead  to  the  formation  of  the  contact  areas  as  seen  in  the  above  figures.  However, 
this  issue  does  not  bear  a  high  significance  for  the  calculation  of  the  SIFs.  For  the 
accuracy  of  the  SIF  values,  the  dominating  variable  is  the  total  contact  load  transferred 
rather  than  the  exact  distribution. 

Figures  3.5  and  3.6  show  the  comparison  of  location  of  contact  loads  by  Hertzian 
analysis  and  FE  contact  analysis  for  both  single-tooth  and  double-tooth  contact  at  initial 
crack.  Gray  squares  in  the  figures  define  the  extreme  points  of  the  Hertzian  contact 
ellipses.  Black  points  are  the  mesh  nodes  where  contact  occurs  in  the  FE  contact  analysis. 
The  location  of  contact  loads  is  almost  the  same  at  crack  #1  and  crack  #2. 


Figure  3.5:  Comparison  of  location  of  contact  loads  by  Hertzian  analysis  and  FE  contact 
analysis  for  double-tooth  contact  at  initial  crack.  Hertzian  loads  on  tooth  1  and  tooth  2 
correspond  to  load  step  3  and  load  step  13,  respectively. 
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Figure  3.6:  Comparison  of  location  of  contact  loads  by  Hertzian  analysis  and  FE  contact 
analysis  for  double-tooth  contact  at  initial  crack.  Hertzian  load  on  tooth  1  corresponds  to 
load  stepll. 

3.3.2  Comparison  of  SIFs  extracted  from  three-dimensional  FEM  contact  analysis  and 
finite  element  analysis  with  Hertz  contact  loads 

Change  in  contact  loads  due  to  cracking  is  expected  to  have  an  influence  on  the 
stress  intensity  factors  and  therefore  on  the  crack  trajectory.  In  order  to  investigate  this 
issue,  SIFs  were  calculated  using  the  results  of  three-dimensional  finite  element  contact 
analysis  and  finite  element  analysis  using  Hertzian  contact  loads.  The  latter  set  of 
analyses  were  carried  out  using  a  two-tooth  pinion  model  where  the  Hertzian  contact 
loads,  corresponding  load  step  3  and  load  step  1 1  of  the  fifteen  discretized  moving  load 
steps  described  in  section  2.1,  were  applied  to  the  model. 

SIFs  were  first  calculated  by  the  displacement  correlation  technique  for  both 
cases;  however,  contact  analysis  results  did  not  compare  well  with  the  Hertzian  contact 
load  results.  A  30  %  difference  in  SIFs  was  observed  between  the  two  set  of  analysis.  We 
used  linear  elements  for  contact  calculations,  which  is  required  for  this  type  of  analysis, 
and  extracted  SIFs  from  these  results  by  the  displacement  correlation  technique.  This 
approach  is  not  accurate  enough  with  the  linear  elements. 

The  equivalent  domain  J-integral  method  is  known  to  give  more  accurate  results, 
and  this  method  was  implemented  in  FRANC3D.  However,  SIFs  obtained  by  linear 
elements  were  still  somewhat  less  accurate  than  the  quadratic  elements  used  for  the  finite 
element  analysis  with  Hertz  contact  loads.  Improved  results  with  the  linear  elements  were 
obtained  as  die  size  of  the  radius  determining  the  domain  of  the  J-integral  calculations 
increased. 

Comparisons  between  Hertz  contact  loads  and  contact  loads  calculated  by  three- 
dimensional  finite  element  analysis  were  given  in  the  previous  section.  In  this  section,  we 
present  the  effects  of  the  difference  in  contact  loads  calculated  by  two  different  methods 
on  the  SIFs.  Figures  3.7  through  3.12  compare  mode  I,  II  and  III  SIFs  extracted  from  the 
results  of  the  analyses  with  Hertz  contact  conditions  and  three-dimensional  finite  element 
contact  analyses. 

SIF  values  were  calculated  for  the  models  mentioned  in  the  previous  section, 
namely,  a  pinion  with  an  initial  crack,  and  with  cracks  corresponding  to  crack  #1  and 
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crack  #2.  Two  loading  conditions  were  investigated,  one  corresponding  to  double  tooth 
contact  and  the  other  corresponding  to  single  tooth  contact,  assumed  to  represent  load 
step  3  and  load  step  1 1,  respectively. 

The  difference  between  SIFs  calculated  by  three-dimensional  finite  element 
contact  analyses  and  analyses  using  Hertz  contact  loads  increased  as  the  crack  length 
increased  due  to  the  enhanced  flexibility  of  the  tooth.  The  change  in  tooth  flexibility  as 
the  crack  grew,  resulted  in  lower  contact  loads  and  therefore  lower  SIFs. 

For  the  gear/pinion  model  with  an  initial  crack  under  double-tooth  contact,  the 
difference  was  not  expected  to  be  that  significant  since  the  initial  crack  was  very  small 
and  did  not  introduce  significant  additional  flexibility  to  the  tooth.  This  behavior  is 
shown  in  Figure  3.7.  At  crack  #1,  as  shown  in  Figure  3.8,  the  change  in  the  SIFs  became 
more  pronounced.  An  even  more  amplified  effect  of  the  change  in  tooth  flexibility  was 
observed  for  crack  #2,  as  plotted  in  Figure  3.9. 


■  KIcontact 
□  KIHertz 

♦  KII_contact 
o  KII__Hertz 

•  KHIcontact 
o  KIII  Hertz 


Figure  3.7:  Comparison  of  mode  I,  II  and  III  SIFs  calculated  by  three-dimensional  finite 
element  contact  analysis  and  finite  element  analysis  with  Hertzian  contact  analysis  for  a 
pinion  with  initial  crack  and  at  a  double  tooth  contact  configuration. 
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normalized  crack  front 


Figure  3.8:  Comparison  of  mode  I,  II  and  III  SIFs  calculated  by  three-dimensional  finite 
element  contact  analysis  and  finite  element  analysis  with  Hertzian  contact  analysis  for  a 
pinion  with  crack  #1  and  at  a  double  tooth  contact  configuration. 


normalized  crack  length 


Figure  3.9:  Comparison  of  mode  I,  II  and  III  SIFs  calculated  by  three-dimensional  finite 
element  contact  analysis  and  finite  element  analysis  with  Hertz  contact  loads  for  a  pinion 
with  crack  #2  and  at  a  double  tooth  contact  configuration. 
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The  analyses  were  repeated  for  an  almost  single  tooth  contact  case  on  all  three  of 
the  crack  configurations.  The  same  trends  were  observed  as  in  the  double  tooth  contact  as 
seen  in  Figures  3.10  through  3.12.  However,  the  difference  between  the  contact  analyses 
and  analyses  with  Hertz  contact  loads  were  higher  for  this  case.  This  might  be  due  to  the 
fact  that  the  total  load  applied  to  the  tooth  was  less  than  the  total  load  applied  in  analysis 
with  Hertzian  contact.  As  described  in  section  3.3.1,  and  shown  in  Table  3.1,  there  was  at 
least  a  very  little  portion  carried  by  the  adjacent  tooth  decreasing  the  amount  of  total  load 
carried  by  the  cracked  tooth. 
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Figure  3.10:  Comparison  of  mode  I,  II  and  III  SIFs  calculated  by  three-dimensional  finite 
element  contact  analysis  and  finite  element  analysis  with  Hertzian  contact  loads  for  a 
pinion  with  initial  crack  and  at  a  single  tooth  contact  configuration. 

Figures  3.13  through  3.15  show  the  variation  of  the  ratio  of  Ku  to  Kj  along  the 
crack  front  for  both  double-tooth  and  single-tooth  contact.  These  graphs  are  plotted  using 
the  results  obtained  by  FE  contact  analysis  for  all  three  models  with  initial  crack,  crack 
#1  and  crack  #2. 
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normalized  crack  front 

Figure  3.11:  Comparison  of  mode  I,  II  and  III  SIFs  calculated  by  three-dimensional  finite 
element  contact  analysis  and  finite  element  analysis  with  Hertzian  contact  loads  for  a 
pinion  with  crack  #1  and  at  a  single  tooth  contact  configuration. 


normalized  crack  length 

Figure  3.12:  Comparison  of  mode  I,  II  and  III  SIFs  calculated  by  three-dimensional  finite 
element  contact  analysis  and  finite  element  analysis  with  Hertzian  contact  loads  for  a 
pinion  with  crack  #2  and  at  a  single  tooth  contact  configuration. 
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normalized  crack  front 

Figure  3.13:  Ratio  of  Kn  to  Kj  along  the  crack  front  for  double-tooth  and  single-tooth 
contact  at  initial  crack. 


Double-tooth  contact 


^  0.6 


Single-tooth  contact 


0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9 


normalized  crack  front 

Figure  3.14:  Ratio  of  Kn  to  K]  along  the  crack  front  for  double-tooth  and  single-tooth 
contact  at  crack  #1 . 
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normalized  crack  front 


Figure  3.15:  Ratio  of  Ku  to  Kj  along  the  crack  front  for  double-tooth  and  single-tooth 
contact  at  crack  #2. 

By  performing  these  contact  analyses,  we  were  able  to  show  that  the  change  in  the 
stiffness  of  the  tooth  due  to  a  growing  crack  has  significant  effects  on  the  SIFs  and 
therefore  on  the  crack  trajectory  and  fatigue  life.  If  LEFM  analyses  of  section  2  are 
repeated  under  these  new  contact  loading  conditions,  the  change  in  the  crack  trajectory 
and  the  fatigue  life  of  the  pinion,  as  a  result  of  taking  into  account  the  change  in  the  tooth 
flexibility,  can  be  determined. 


4.  Summary 

This  report  summarized  new  results  for  predicting  crack  trajectory  and  fatigue  life 
for  a  spiral  bevel  pinion  using  the  Finite  Element  Method  (FEM).  The  predictions 
presented  in  this  report  are  based  on  LEFM  theories  combined  with  the  FEM, 
incorporating  plasticity  induced  fatigue  crack  closure  and  moving  loads. 

Previously  in  this  project,  three-dimensional  computational  simulation  of  crack 
growth  in  spiral  bevel  gears  was  done  using  BEM  [1].  As  a  result  of  these  analyses,  a 
method  for  predicting  three-dimensional,  non-proportional,  fatigue  crack  growth 
incorporating  moving  loads  was  developed.  The  BEM  studies  showed  the  need  for 
improving  the  efficiency  and  the  accuracy  of  the  computations  by  employing  the  state-of- 
the-art  computational  capabilities. 

In  the  current  work,  three-dimensional  finite  element  modeling  of  the  spiral  bevel 
pinion  was  done  using  OSM/FRANC3D.  A  three-dimensional  FEM  mesh  of  the  pinion 
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model  was  created  by  a  3D  mesher.  Moving  contact  loads  on  the  tooth  surface  were 
calculated  by  interpolation  using  the  shape  functions  on  the  surfaces  of  the  loaded 
elements.  The  analyses  were  carried  out  using  a  parallel  FEM  solver,  which  calculates 
SIFs  using  equivalent  domain  J-integral  method.  Fatigue  life  predictions  were  made 
based  on  a  modified  Paris  model  incorporating  crack  closure. 

In  this  report,  we  show  that  we  can  simulate  fatigue  crack  growth  in  a  spiral  bevel 
gear  more  accurately  and  efficiently  using  the  FEM  along  with  a  better  representation  of 
moving  loads.  Another  very  significant  improvement  was  the  decrease  in  solving  time  of 
the  problem  by  employing  parallel  PC-clusters.  This  reduces  the  computation  time  to  a 
few  minutes  while  substantially  increasing  the  resolution  of  the  solution. 

To  obtain  a  more  detailed  understanding  of  the  contact  between  a  cracked  pinion 
tooth  in  mesh  with  an  uncracked  gear  tooth,  three-dimensional  contact  analyses  were 
performed  on  a  spiral  bevel  gear  set  incorporating  a  crack.  The  goal  in  carrying  out  these 
analyses  was  to  capture  the  redistribution  of  contact  loads  due  to  crack  growth.  Results  of 
these  analyses  showed  the  expected  trend  of  decreasing  tooth  loads  carried  by  the  cracked 
tooth  with  increasing  crack  length.  We  also  showed  that  this  decrease  in  contact  loads 
had  an  impact  on  the  SIF  values  and  therefore  would  also  affect  the  crack  trajectory  and 
fatigue  life  predictions. 

Predicting  crack  trajectories  more  efficiently  and  accurately  is  very  important  for 
gear  designers.  Development  of  practical  and  accurate  numerical  tools  for  evaluating  the 
gear  performance  and  life  enables  for  new  and  better  gear  designs.  The  work  presented 
here  provides  such  a  tool  capable  of  predicting  crack  trajectory  and  fatigue  life  in  a  spiral 
bevel  gear. 
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Appendix:  Abbreviated  sample  ABAQUS  input  file  for  contact  analysis 


* HEADING 
FRANC3D  Model 

**  FRANC3D  Version  2.0  Generated  Input  File 


**  NUMBER  OF  NODES:  13029 
**  NUMBER  OF  ELEMENTS:  58712 

*  * 

**  NODE  DEFINITION 

♦NODE 

1, 

3.815870,  0.262049, 

1.500000 

2, 

3.585600,  0.038081, 

1.500000 

3, 

3.824640,  0.040620, 

1.500000 

13028, 

5.109482,  0.491277, 

0.734872 

13029, 

5.110899,  0.478949, 

0.730531 

**  END  NODE  DEFINITION 
** 

**  ELEMENT  DEFINITION 

♦ELEMENT,  TYPE=C3D4,  ELSET=PID1 

1,1904,1855,2100,1903 

* ELEMENT,  TYPE=C3D4,  ELSET=PID1 

2,1904,2100,1953,1952 

♦ELEMENT,  TYPE=C3D4,  ELSET=PID1 

3,125,107,106,846 


♦ELEMENT,  TYPE=C3D4,  ELSET=PID2 

58711.7717.12943.7703.13029 
♦ELEMENT,  TYPE=C3D4,  ELSET=PID2 

58712. 6078.7717.7703.13029 
**  END  ELEMNT  DEFINITION 
** 

*SOLID  SECTION,  ELSET=PID1,  MATERI AL=MAT  1 
*  MATERIAL,  NAME=MAT1 
^ELASTIC,  TYPE=ISO 
3e+07 ,  0.3 
^DENSITY 
7800 

*S0LID  SECTION,  ELSET=PID2,  MATERI AL=MAT 2 
* MATERIAL,  NAME=MAT2 
^ELASTIC,  TYPE=ISO 
3e+07,  0.3 
^DENSITY 
7800 
** 

**Definition  of  contact  for  tooth  1 


*ELSET ,  ELSET=MASTER_1 
26454, 

26457, 

26458, 
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57531 , 

58366, 

58367, 

*ELSET,  ELSET=SLAVE_1 

1, 

2, 

5, 


22208, 

24437, 

24438, 

* SURFACE  DEFINITION,  NAME=MASTER_1 

26454. 54 

26457. 54 

26458. 54 


57531. 51 

58366. 51 

58367. 51 

* SURFACE  DEFINITION,  NAME=SLAVE_1 

1,  S2 

2,  S4 
5,  S4 


22208. 51 

24437. 51 

24438. 51 

* SURFACE  INTERACTION,  NAME=NOFRIC 
* FRICTION 
0 

* CONTACT  PAIR,  INTERACTION-NOFRIC 
MASTER_1,  SLAVE__1 
** 

**  Definition  of  contact  for  tooth  2 
** 

*ELSET ,  ELSET=MASTER__2 
26455, 

26456, 

26459, 


57320, 

57753, 

57754, 

*ELSET ,  ELSET=SLAVE_2 

3, 

4, 

9, 
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24532, 

24540, 

24548, 

* SURFACE  DEFINITION,  NAME=MASTER_2 

26455. 54 

26456. 54 

26459. 54 


57320. 51 

57753. 51 

57754. 51 

* SURFACE  DEFINITION,  NAME=SLAVE_2 

3,  S2 

4 ,  S4 
9,  S4 


24532. 51 

24540. 51 

24548. 51 

*SURFACE  INTERACTION,  NAME=NOFRICl 
*FRICTION 
0 

*CONTACT  PAIR,  INTERACTION=NOFRICl 

MASTER_2,  SLAVE_2 
** 

**  Application  of  displacement  to  ensure  contact 
**  between  gear  and  pinion 


*STEP,  NLGEOM,  INC=100 
**FRANC3D  load  case  1 
*STATIC 


1.0,  1.0 

* BOUNDARY 

2209.3.3,  0 

2211.3.3,  0 

2213.3.3,  0 

2215.3.3,  0 

2217.3.3,  0 

2219.3.3,  0 

2221.3.3,  0 

2223.3.3,  0 

2232.2.2,  -0.01 

2247.3.3,  0 


33,1,1,0 

33,2,2,0 

33,3,3,0 

2322,1,1,0 
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2322 ,2,2,0 
2322 ,3,3,0 
2415/ 1,1,0 
2415,2,2/0 
2415,3,3,0 
*END  STEP 
** 

**  Application  of  full  design  torque 

*STEP,  NLGEOM,  INC=100 
♦♦FRANC3D  load  case  1 
♦STATIC 
1.0,  1.0 

♦BOUNDARY,  OP=NEW 
.2209,3,3,  0 

2211.3.3,  0 

2213.3.3,  0 

2215.3.3,  0 

2217.3.3,  0 

2219.3.3,  0 

2221.3.3,  0 

2223.3.3,  0 

2247.3.3,  0 


33,1,1,0 

33,2,2,0 

33,3,3,0 

2322,1,1,0 

2322,2,2, 0 

2322,3,3,0 

2415,1,1,0 

2415,2,2,0 

2415,3,3,0 

♦CLOAD 

2799.2, -1000 
♦END  STEP 

♦STEP,  NLGEOM,  INC=100 
♦♦FRANC3D  load  case  1 
♦STATIC 
1.0,  1.0 
♦CLOAD 

2799.2, -3230 
♦FILE  FORMAT,  ASCII 
♦NODE  FILE,  FREQUENCY=1 0 0 
U 

♦EL  FILE,  FREQUENCY=1 0 0 ,  POSITION=NODES 
S,  E,  PE 

♦RESTART, WRITE,  FREQUENCY=100 
♦NODE  PRINT,  FREQUENCY=1 0 0 
U 

♦EL  PRINT,  FREQUENCY=100,  POSITION=NODES 

♦CONTACT  PRINT, FREQUENCY=1 00 

CFN, CARE A, XN 

♦END  STEP 

♦♦END  OF  INPUT 
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